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Likelihood ratio tests are performed for the hypothesis that charged-particle multiplicities mea- 
sured in Au-Au and Cu-Cu collisions at tJsnn = 62.4 and 200 GeV are distributed according to 
the negative binomial form. Results indicate that the hypothesis should be rejected in the all cases 
of PHENIX-RHIC measurements. Possible explanations of that and of the disagreement with the 
least-squares fitting method are given. 
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I. INTRODUCTION 

The analysis of charged hadron multiplicities in Au-Au and Cu-Cu collisions at ^/snn = 62.4 and 200 GeV was 
done by the PHENIX Collaboration in [l| . It was also claimed there that these multiplicities are distributed according 
to the negative binomial form. The UA5 Collaboration noticed for the first time that charged-particle multiplicity 
distributions measured in high energy proton- (anti)proton collisions in limited intervals of pseudo-rapidity have this 
form 0. 

The Negative Binomial Distribution (NBD) is defined as 

P(n; P ,k) = Kk + D(k + 2) (k + n-1) k 

where n = 0, 1, 2, < p < 1 and k is a positive real number. In the application to high energy physics n has the 
meaning of the number of charged particles detected in an event. The expected value n and variance V(n) [l~o| are 
expressed as: 

n = , Vin) = . (2) 

p p z 

Multiplicity fluctuations are expressed in terms of the scaled variance: 



(N 2 ch ) - (N ch )* V(n) 
(N ch ) n 



(3) 



where N c h is the charged particle multiplicity and the last equality is valid only for the whole population (the set of 
all possible outcomes if the experiment is repeated infinitely many times), assuming that the hypothesis about the 
NBD is true. 

In application to the high energy physics, the parameters fc,n instead of k,p are used usually and 



which is the scaled variance, Eq. ([3]). But because the centrality bins have the nonzero width, fluctuations defined by 
Eq. ([3]) also include a non-dynamical component. This component is the result of the fluctuations of the geometry of the 
collisions within a given centrality bin. The geometrical fluctuations were evaluated by the PHENIX Collaboration 
in It turned out that those fluctuations can be expressed by a correction factor, f geo , which is independent 
of centrality but varies with the collision type. Then the pure scaled variance now representing only dynamical 
fluctuations, i.e. after subtraction of the geometrical component, can be calculated from the following equation [l|: 



l^dyn — 1 — fgeo 1 — 1) 

Also parameter k changes to k& yn accordingly 



(5) 
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^dyn — fg e ° ' ^ 1 ■ (6) 

In this analysis the hypothesis that the charged-particle multiplicities measured in ultra-relativistic heavy-ion colli- 
sions are distributed according to the NBD is verified with the use of the maximum likelihood method (ML) and the 
likelihood ratio test. More details of this approach can be found in Refs. 

There are two crucial reasons for this approach: 

1. The fitted quantity is a probability distribution function (p.d.f.), so the most natural way is to use the ML 
method, where the likelihood function is constructed directly from the tested p.d.f.. In fact, what is fitted are 
parameters of the distribution. The fitted values are the estimators of these parameters. It is well-known in 
mathematical statistics that an ML estimator is consistent, asymptotically unbiased and efficient [!, 0, Hj] . But 
even more important is that because of Wilks's theorem (see Appendix B) one can easily define a statistic, the 
distribution of which converges to a \ 2 distribution as the number of measurements goes to infinity. Thus for the 
large sample the goodness-of-fit can be expressed as a p- value computed with the corresponding \ 2 distribution. 

2. The most commonly used method, the least-squares (LS) method (called also the x 2 minimization), has the 
disadvantage of providing only the qualitative measure of the significance of the fit, in general. Only if observables 
are represented by Gaussian random variables with known variances, the conclusion about the goodness-of-fit 
equivalent to that mentioned in the point 1 can be derived Q. 

It is worth noting that the ML method with binned data and Poisson fluctuations within a bin was already applied 
to fitting multiplicity distributions to the NBD but at much lower energies (E-802 Collaboration 0). 



II. LIKELIHOOD RATIO TEST 



The number of charged particles N c h is assumed to be a random variable with the p.d.f. given by Eq. (fTJ). Each 
event is treated as an independent observation of and a set of a given class of events is a sample. For N events 
in the class there are N measurements of N c h, say X = {X±, X2, X^}. Some of these measurements can be equal, 
i.e. Xi = Xj for i ^ j can happen. The whole population consists of all possible events with the measurements of 0, 
1, 2,... charged particles and by definition is infinite (Tlj . 

Let divide the sample into m bins characterized by Yi - the number of measured charged particles [l2| and rii - the 
number of entries in the ith bin, N — ^2 i=1 n.i (details of the theoretical framework of this Section can be found in 
Refs. 0-0] )■ Then the expectation value of the number of events in the ith bin can be written as 

Vi(vtot,P,k) = v tot ■ P(Y l ;p,k) , (7) 

where vtot is the expected number of all events in the sample, vtot — ^2i—i v i- This is because one can treat the 
number of events in the sample N also as a random variable with its own distribution - Poisson one. Generally, the 
whole histogram can be treated as one measurement of m-dimensional random vector n = (m, ...,n TO ) which has a 
multinomial distribution, so the joint p.d.f. for the measurement of N and n can be converted to the form (3,0): 

m v ni 

/(n; 1/1, v m ) = J[ J- exp (-^) . (8) 
Since now /(n; v m ) is the p.d.f. for one measurement, / is also the likelihood function 

L(n I v x , v m ) = /(n; v x , v m ) . (9) 
With the use of Eq. ([7|) the corresponding likelihood function can be written as 



L(n I v toU p,k) = L(n | vi(v totl p,k), ...,u m {v tot ,p,k)) 
Then the likelihood ratio is defined as 



(10) 
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A _ L(n | i> to t,P,k) _ L(n \ v t ot,P,k) 
L{n | vx,...,v m ) L(n | m,...,n m ) ' 

where Vtot, P an d are the ML estimates of vtot, P and fe with the likelihood function given by Eq. (|10p and i/j = rij, 
i = 1, 2, ...m are the ML estimates of i/j treated as free parameters. Note that since the denominator in Eq. (fTTj) does 
not depend on parameters, the log-ratio defined as 

i w mi L(n\vtot,P,k) 

m\{u t ot,P, k) = In ■ 



L(n | m, ...,n n 



m 



Tij In V Vi — n., 

»=i 



m 

n 



AT - ^ n< In— , (12) 



where Vi are expressed by Eq. {?]), can be used to find the ML estimates of v to t, P and k. The values z> tot , p and A; for 
which \{v to t,p 1 k) has its maximum are the maximum likelihood estimates of parameters vtot, P and k. Further, the 
statistic given by 



i=l 



X 2 = -2ln\(v tot ,p, k) = 2^ ( Uj ln^ + A; - rn ) (13) 



approaches the \ 2 distribution asymptotically, i.e. as the number of measurements, here the number of events N , 
goes to infinity (the consequence of the Wilks's theorem, see Appendix B). The values z)j are the estimates of Vi given 
by 

h = f>tofP(Yi;p,k) (14) 
and if one assumes that u to t does not depend on p and k then vt ot — N. For such a case 

m rn 



and Eq. (|T3|) becomes 



i=l i=l 



m 

X 2 (p,fc) = -2 In A = 2 V ri 4 ln— . (16) 



Also then one can just put v tot — N and Eq. (1121) can be rewritten as 



In A(p, k) 

m m 

= N -hxN - ^ n-ilnni + X n l \nP{Y i ]p, k) 

i=l i=l 
m m 

i=l i=l 

m m 

= -ny, pr^pr+Nj2 pr^p(Yi-,p,k), (n) 

i=l i=l 
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where Pf x = rii/N. Note that the maximum of In A is the minimum of \ 2 = ~2 m A, so from Eqs. ([To]) and p7|) one 
arrives at 



x ^ n = -2Arf;^ln^glM. (18) 
It can be proven that one of the necessary conditions for the existence of the maximum is (see Appendix A for details) : 



fi = (N ch ) , (19) 

i.e. the distribution average has to be equal to the experimental average. This is very good because (N c h) is what is 
called in statistics a sample mean. The sample mean is an estimator for the expectation value of the random variable, 
which is consistent and unbiased [3]. In other words the ML estimator of n is (N c h) (n — (N c h)). 



III. RESULTS AND DISCUSSION 



The method described in Section [TTI requires that all bins in a given data set have the width equal to 1, so as the 
experimental probability Pf x to measure a signal in the ith bin was equivalent to the probability of the measurement 
of (i — 1) charged particles (the first bin is the bin of charged particles detected). This is fulfilled for all bins of the 
considered data sets 

Since the test statistic —2 In A has a \ 2 distribution approximately in the large sample limit, it can be used as a 
test of the goodness-of-fit. The result of the test is given by the so-called p- value which is the probability of obtaining 
the value of the statistic, Eq. (TT5T) . equal to or greater then the value just obtained for the present data set, when 
repeating the whole experiment many times: 



P = P(X 2 > Xminl n d) = \ fi*] n d)dz , (20) 

J x 2 min 

where f(z] rid) is the \ 2 p.d.f. and nd the number of degrees of freedom, nd — m — 2 here. 

The results of the analysis are presented in Tables Hl lVIIII and illustrated with Figs. [TJSJ In fact the whole analysis 
was done for the two kinds of histograms: (i) bins with the number of entries rii < 5 excluded, Tables HI IIII1 IVland 
IVIII (ii) bins with the number of entries rij < 60 excluded, Tables [TTI IIV1 IVTl and IVIII1 In practice this corresponds 
to cutting off less (i) or more (ii) the tails of the full measured histograms. The tails breaks the visual agreement 
between the data and the NBD, cf. Figs. Q] and [5J The condition that only bins with m > 5 are taken into account 
is the minimal condition imposed on a histogram to do any statistical inference without Monte Carlo simulations Q . 
The condition (ii) corresponds to the choice made originally by the PHENIX Collaboration in their analysis [l|. It 
has turned out that the results of this analysis are qualitatively the same for both choices. 

As one can see, the hypothesis in question should be rejected in all considered cases. But it was claimed that 
charged-particle multiplicities measured in Au-Au and Cu-Cu collisions at ^/snn — 62.4 and 200 GeV are distributed 
according to the NBD However that conclusion was the result of the application of the LS method. Therefore 
it seems to be reasonable to check what are the values of the LS test statistic at the ML estimators listed in the 
third and fourth columns of Tables Hl lVIIII For the sample described in Sect. [ITJone can define the LS test statistic 
(commonly called the \ 2 function) as: 



*UM) -f wr-J-cM . (21) 

i—l 1 

where erri is the uncertainty of the ith measurement. Here this function is not minimized with respect to n and 
k (or p and k) as in the LS method but is calculated at ML estimates of n and k. This is allowed in statistics and 
is equivalent to test a single point in the parameter space. Then the tested point might not be the best estimate of 
the true value but the hypothesis in question becomes the hypothesis only about a particular distribution (a simple 
hypothesis). One can see from the ninth column of Tables Hl lVIIII that XLsl n d va hies are significant for almost all 
centrality classes, what agrees with the results of Ref. 1]. But this contradicts the results of the likelihood ratio test, 
which are expressed by x 2 / n d & n d p- values listed in the seventh and eight columns of Tables Hl lVIIII 
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TABLE I: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Au-Au collisions at ^/snn = 
200 GeV, f g eo = 0.37 ± 0.027 1]. Fitting ranges are limited to the bins with m > 5, where rii is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 






x 2 

( n d) 


p- value 

[%] 


xls/ n d with errors: 
quadrature statistical 
sum only 


0-5 


653145 


270.0 
±2.5 


61.85 
±0.01 


1.37 -10~ 3 
±0.10 -10" 3 


1.08 
±0.01 


23.73 
1756.0 
(74) 





0.98 9.8 


5-10 


657944 


163.4 
±1.2 


53.91 
±0.01 


2.26 -10" 3 
±0.17 -10 -3 


1.12 
±0.01 


9.12 
592.7 
(65) 





0.69 6.9 


10-15 


658739 


112.5 
±0.7 


46.50 
±0.01 


3.29 -10~ 3 
±0.24 -lO -3 


1.15 
±0.01 


11.5 
795.5 
(69) 





0.66 6.6 


15-20 


659607 


85.1 
±0.5 


39.72 
±0.01 


4.35 -10" 3 
±0.32 -10 -3 


1.17 
±0.01 


8.9 

585.8 
(66) 





0.52 5.2 


20-25 


658785 


67.6 
±0.4 


33.56 
±0.01 


5.48 -10~ 3 
±0.40 -10 -3 


1.18 
±0.01 


13.5 
848.8 
(63) 





0.46 4.6 


25-30 


659632 


56.7 
±0.3 


28.01 
±0.01 


6.52 -10" 3 
±0.48 -10" 3 


1.18 
±0.01 


10.9 
640.6 
(59) 





0.37 3.7 


30-35 


659303 


47.4 
±0.3 


23.02 
±0.01 


7.81 -10" 3 
±0.57 -10~ 3 


1.18 
±0.01 


7.9 
429.9 

(54) 





0.31 3.1 


35-40 


661174 


40.5 
±0.2 


18.64 
±0.01 


9.13 -10~ 3 
±0.67 -10" 3 


1.17 
±0.01 


8.5 
389.7 
(46) 





0.37 3.7 


40-45 


661599 


34.0 
±0.2 


14.84 
±0.01 


1.09 -10" 2 
±0.80 -lO" 3 


1.16 
±0.01 


7.3 
301.0 

(41) 





0.35 3.5 


45-50 


661765 


27.3 
±0.2 


11.57 
±0.005 


1.35 -10" 2 
±0.99 -lO" 3 


1.16 
±0.01 


10.5 
390.2 
(37) 





0.92 9.2 


50-55 


662114 


21.3 
±0.1 


8.82 
±0.004 


1.74 -lO -2 
±0.13 -10" 2 


1.15 
±0.01 


38.8 
1436.4 





12.06 120.6 
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TABLE II: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Au-Au collisions at tJsnn = 
200 GeV, fgeo = 0.37 ± 0.027 Q. Fitting ranges are limited to the bins with ru > 60, where ru is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 






x 2 

( n d) 


p- value 

[%] 


xls/ n d with errors: 
quadrature statistical 
sum only 


0-5 


652579 


289.0 
±2.9 


61.86 
±0.01 


1.28 -10~ 3 
±0.94 -10 -4 


1.08 
±0.01 


20.0 
1160.2 

(58) 





0.57 5.7 


5-10 


657571 


168.1 
±1.2 


53.91 
±0.01 


2.20 -10 -3 
±0.16 -10~ 3 


1.12 
±0.01 


20.56 
1151.6 

(56) 





0.61 6.1 


10-15 


658258 


116.4 
±0.7 


46.50 
±0.01 


3.18 -10~ 3 
±0.23 -10 -3 


1.15 
±0.01 


18.4 
991.7 

(54) 





0.53 5.3 


15-20 


659302 


86.9 
±0.5 


39.72 
±0.01 


4.26 -10 -3 
±0.31 -10 -3 


1.17 
±0.01 


12.6 
667.5 
(53) 





0.43 4.3 


20-25 


658461 


69.1 
±0.4 


33.56 
±0.01 


5.36 -lO -3 
±0.39 -10~ 3 


1.18 
±0.01 


12.3 
604.7 

(49) 





0.34 3.4 


25-30 


659337 


57.9 
±0.3 


28.0 
±0.01 


6.39 -10" 3 
±0.47 -10" 3 


1.18 
±0.01 


10.4 
469.1 

(45) 





0.28 2.8 


30-35 


659021 


48.3 
±0.3 


23.02 
±0.01 


7.66 -10" 3 
±0.56 -10" 3 


1.18 
±0.01 


8.6 
351.02 
(41) 





0.16 1.6 


35-40 


660937 


41.3 
±0.2 


18.64 
±0.01 


8.96 -10~ 3 
±0.66 -10 -3 


1.17 
±0.01 


7.6 

280.3 
(37) 





0.19 1.9 


40-45 


661422 


34.6 
±0.2 


14.84 
±0.01 


1.07 -10" 2 
±0.78 -lO" 3 


1.16 
±0.01 


7.9 
260.3 
(33) 





0.21 2.1 


45-50 


661577 


27.9 
±0.2 


11.56 
±0.005 


1.33 -10~ 2 
±0.97 -lO" 3 


1.15 
±0.01 


10.0 
279.9 
(28) 





0.23 2.3 


50-55 


661877 


21.9 
±0.1 


8.81 
±0.004 


1.69 -lO" 2 
±0.12 -10" 2 


1.15 
±0.01 


40.0 
959.2 





0.30 3.0 
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FIG. 1: Uncorrected multiplicity distributions of charged hadrons for 200 GeV Au-Au collisions [l( within ranges limited to 
the bins with nj > 5. The lines are fits to the NBD. The data are scaled by the amounts in the legend. Errors represent the 
statistical and systematic errors added in quadrature. 
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FIG. 2: Uncorrected multiplicity distributions of charged hadrons for 200 GeV Au-Au collisions [l| within ranges limited to 
the bins with m > 60. The lines are fits to the NBD. The data are scaled by the amounts in the legend. Errors represent the 
statistical and systematic errors added in quadrature. 



However, the ML method and the likelihood ratio test have much more solid mathematical foundations as it has 
been explained in Sect. |U Thus the hypothesis about the NBD should be rejected on the basis of obtained values of 
X 2 /n<i and p- values. 

The crucial question is now why the conclusions from \ 2 an d Xls * es ^ statistics are entirely opposite for PHENIX 
measurements? The main difference between both statistics is that \ 2 depends explicitly on the number of events 
but Xls does not. On opposite, \ 2 does not depend on the actual errors but x\s does. In fact, x 2 statistic implicitly 
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TABLE III: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Au-Au collisions at ^/snn = 
62.4 GeV, f geo = 0.33 ± 0.031 [lj. Fitting ranges are limited to the bins with n; > 5, where rii is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 


1 / kdyn 




X 2 /n d 

x 2 

(rid) 


p- value 

[%] 


XLs/n>d with errors: 
quadrature statistical 
sum only 


0-5 


607155 


225.2 
±2.5 


44.67 
±0.01 


1.47 -10~ 3 
±0.14 -10~ 3 


1.07 
±0.01 


2.37 
139.6 

(59) 


1.7-10" 8 


0.18 1.8 


5-10 


752392 


142.3 
±1.1 


37.96 
±0.01 


2.32 TO -3 
±0.22 -10 -3 


1.09 
±0.01 


2.44 
131.9 

(54) 


1.9-10" 8 


0.11 1.1 


10-15 


752837 


115.2 
±0.9 


31.53 
±0.01 


2.87 -10" 3 
±0.27 -10" 3 


1.09 
±0.01 


2.06 
107.1 

(52) 


1.1-10" 5 


0.13 1.3 


15-20 


752553 


88.0 
±0.6 


26.07 
±0.01 


3.75 TO -3 
±0.35 TO" 3 


1.10 
±0.01 


1.86 
87.3 
(47) 


3.2-10" 4 


0.13 1.3 


20-25 


752296 


68.5 
±0.5 


21.35 
±0.01 


4.82 -10~ 3 
±0.45 -10 -3 


1.10 
±0.01 


2.63 
113.2 

(43) 


3.1-10" 8 


0.21 2.1 


25-30 


752183 


53.2 
±0.4 


17.30 
±0.01 


6.21 -10~ 3 
±0.59 TO -3 


1.11 
±0.01 


2.75 
107.3 

(39) 


2.7-10~ 8 


0.23 2.3 


30-35 


751375 


40.1 
±0.3 


13.84 
±0.005 


8.22 TO -3 
±0.77 -10~ 3 


1.11 
±0.01 


2.97 
103.9 

(35) 


9.6-10" 9 


0.25 2.5 


35-40 


751661 


31.7 
±0.2 


10.89 
±0.004 


1.04 TO" 2 
±0.98-10 


1.11 
±0.01 


6.72 
194.9 
(29) 





0.16 1.6 


40-45 


750884 


25.1 
±0.2 


8.42 
±0.004 


1.31 TO -2 
±0.12 -10~ 2 


1.11 

±0.01 


37.5 
937.4 
(25) 





40.36 403.6 


45-50 


751421 


21.8 
±0.2 


6.41 
±0.003 


1.51 -10~ 2 
±0.14 -10~ 2 


1.10 
±0.01 


209.0 
4806.8 
(23) 





285.9 2859.5 



assumes errors of the type ,/n7, what is the straightforward result of the form of the likelihood function, Eqs. ([5]) and 
©, namely the product of Poisson distributions. In the PHENIX analysis [l| errors erri in Eq. (|2Tj) are represented 
by the quadrature sum of the statistical and systematic components. And the systematic errors were estimated as 
3.0 times the statistical errors Q. This causes that errf = 10 • <J 2 stat i , where <J s tat,i is the statistical error of the zth 
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FIG. 3: Uncorrected multiplicity distributions of charged hadrons for 62.4 GeV Au-Au collisions [l| within ranges limited to 
the bins with m > 5. The lines are fits to the NBD. The data are scaled by the amounts in the legend. Errors represent the 
statistical and systematic errors added in quadrature. 




FIG. 4: Uncorrected multiplicity distributions of charged hadrons for 62.4 GeV Au-Au collisions [lj] within ranges limited to 
the bins with m > 60. The lines are fits to the NBD. The data are scaled by the amounts in the legend. Errors represent the 
statistical and systematic errors added in quadrature. 



measurement. Hence the values of x\sl nd m ^ as ^ columns of Tables [Ti l Villi are 10 times greater than in preceding 
columns. 

In principle, the accuracy with which experimental distributions approximate the NBD should increase with the 
sample size because if the hypothesis is true the postulated form of distribution is exact for the whole population. 
So with the growing number of events, the experimental distribution should be closer to the postulated one. This is 
also seen in the form of Xmim Eq. ([T5)l . where the linear dependence on N is explicit. To keep Xmin at least constant 
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TABLE IV: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Au-Au collisions at tJsnn = 
62.4 GeV, f geo = 0.33 ± 0.031 [Ij. Fitting ranges are limited to the bins with m > 60, where m is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 






x 2 

( n d) 


p- value 

[%] 


xls/ n d with errors: 
quadrature statistical 
sum only 


0-5 


607075 


227.9 
±2.5 


44.67 
±0.01 


1.45 TO -3 
±0.14 -10" 3 


1.06 
±0.01 


5.55 
294.3 
(53) 





0.19 1.9 


5-10 


752263 


143.9 
±1.1 


37.96 
±0.01 


2.29 TO -3 
±0.22 -10 -3 


1.09 
±0.01 


7.80 
382.4 
(49) 





0.12 1.2 


10-15 


752739 


116.2 
±0.9 


31.53 
±0.01 


2.84 -10~ 3 
±0.27 -10" 3 


1.09 
±0.01 


5.67 
260.8 
(46) 





0.13 1.3 


15-20 


752492 


88.5 
±0.6 


26.07 
±0.01 


3.73 TO -3 
±0.35 TO" 3 


1.10 
±0.01 


5.97 
250.9 
(42) 





0.11 1.1 


20-25 


752182 


69.2 
±0.5 


21.35 
±0.01 


4.77 -10~ 3 
±0.45 -10 -3 


1.10 
±0.01 


10.2 
377.2 

(37) 





0.22 2.2 


25-30 


752095 


53.6 
±0.4 


17.30 
±0.01 


6.16 -10~ 3 
±0.58 -10" 3 


1.11 
±0.01 


8.2 
279.2 
(34) 





0.23 2.3 


30-35 


751324 


40.3 
±0.3 


13.84 
±0.005 


8.19 TO -3 
±0.77 -10~ 3 


1.11 

±0.01 


7.40 
229.3 
(31) 





0.26 2.6 


35-40 


751639 


31.8 
±0.2 


10.89 
±0.004 


1.04 -10" 2 
±0.98 -10" 3 


1.11 

±0.01 


9.43 
254.7 
(27) 





0.15 1.5 


40-45 


750852 


25.2 
±0.2 


8.42 
±0.004 


1.31 TO -2 
±0.12 -10 -2 


1.11 

±0.01 


50.7 
1166.3 

(23) 





0.22 2.2 


45-50 


751348 


22.0 
±0.2 


6.41 
±0.003 


1.50 TO -2 
±0.14 -10 -2 


1.10 
±0.01 


259.8 
4936.4 





343.1 3431.1 



(19) 



when N (the sample size) is growing the relative differences between P(Yi) and Pf x have to decrease. The results of 
the likelihood ratio test have revealed that the discrepancies between the measured and postulated distributions are 
too big for these sample sizes to accept the NBD hypothesis. 

Another surprising point is the application of the LS test statistic with the ML estimators of the parameters to the 
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TABLE V: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Cu-Cu collisions at tJsnn = 
200 GeV, f g eo = 0.40 ± 0.047 [1(. Fitting ranges are limited to the bins with m > 5, where rii is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 


1 / fcdyn 




X 2 /n d 

2 

X 

(n-d) 


p- value 

[%] 


xls/ n d with errors: 
quadrature statistical 
sum only 


0-5 


368510 


59.6 
±0.6 


19.80 
±0.01 


6.72 -10" 3 
±0.79 -lO -3 


1.13 
±0.02 


94.8 
3887.0 
(41) 





2.1 21.2 


5-10 


369206 


49.6 
±0.5 


16.74 
±0.01 


8.06 -10 -3 
±0.95 -10~ 3 


1.13 
±0.02 


16.5 
628.5 
(38) 





0.66 6.6 


10-15 


369945 


41.5 
±0.4 


14.05 
±0.01 


9.64 -10" 3 
±0.11 -10" 2 


1.14 
±0.02 


6.8 
225.5 
(33) 





0.38 3.8 


15-20 


370066 


34.5 
±0.3 


11.78 
±0.01 


1.16 -10" 2 
±0.14 -10 -2 


1.14 
±0.02 


3.0 
92.0 
(31) 


5.8-10" 8 


0.24 2.4 


20-25 


371877 


29.2 
±0.3 


9.81 
±0.01 


1.37 -10" 2 
±0.16 -10 -2 


1.13 
±0.02 


6.6 
186.0 

(28) 





3.4 33.6 


25-30 


368876 


24.9 
±0.2 


8.14 
±0.01 


1.60 -lO" 2 
±0.19 -10" 2 


1.13 
±0.02 


19.3 
502.4 
(26) 





11.5 114.9 


30-35 


368072 


21.9 
±0.2 


6.72 
±0.005 


1.83 -10" 2 
±0.22 -10~ 2 


1.12 
±0.01 


65.6 
1704.8 
(26) 





42.3 422.5 



data samples. The results are presented in the last two columns of Tables Hl l Villi For the choice (ii) which corresponds 
to PHENIX analysis [l|, Tables [TTl IIV1 IVT1 and IVIII1 the XLs/ Ud va l ues obtained here are lower than corresponding 
ones in Ref. pj. Values of the parameters k,n are also different from those in Ref. [l|, what has resulted in slightly 
different (1 — 3% lower) values of the scaled variance ojd yn , see Figs. [7]and[H To make the comparison easier also 
values of kjyn are presented in the fifth column of Tables Hl l Villi Generally, n is greater but the difference does not 
exceed 10% and decreases with the centrality. fc^ n is smaller, especially for case (ii) and the difference also decreases 
with the centrality; from about 20 — 30% for the least central classes to about 5 — 10% for the most central ones. 

IV. CONCLUSIONS 

The main conclusion is that the hypothesis of the NBD of charged-particle multiplicities measured by the PHENIX 
Collaboration in Au-Au and Cu-Cu collisions at y / ijvjv = 62.4 and 200 GeV should be rejected for all centrality 
classes. This is the result of likelihood ratio tests performed for the corresponding data samples. The significant 
systematic errors are the reasons for acceptable values of the least squares test statistic, which is easily seen from the 
last two columns of Tables Hl lVIIII If the systematic errors are not taken into account, then xis/ n d values are greater 
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FIG. 5: Uncorrected multiplicity distributions of charged hadrons for 200 GeV Cu-Cu collisions [lj within ranges limited to 
the bins with m > 5 (left) and m > 60 (right). The lines are fits to the NBD. The data are scaled by the amounts in the 
legend. Errors represent the statistical and systematic errors added in quadrature. 




10 20 30 40 5 10 15 20 25 30 

Nch N ch 

FIG. 6: Uncorrected multiplicity distributions of charged hadrons for 62.4 GeV Cu-Cu collisions [lj within ranges limited to 
the bins with m > 5 (left) and n; > 60 (right). The lines are fits to the NBD. The data are scaled by the amounts in the 
legend. Errors represent the statistical and systematic errors added in quadrature. 



than 1 in all considered cases and substantially greater than 1 in most of them. 

The second, more general conclusion is that the scale of errors should be somehow related to the sample size to make 
the statistical inference trustworthy in the case of the LS method. This is because too big or too small errors cause the 
false inference in this case. Hence errors of the type ^Jnl could be " a frame of reference" . Note that both types of errors 
are such in all considered PHENIX cases; for rn statistical errors are equal to yjnl exactly and systematic errors are 
3.0 times the statistical ones [H,|l], but the last factor causes that the total error is y/l0-(J s tat,i ~ 3.0-a sta t,i — 3.0 -y 7 ^. 
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TABLE VI: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Cu-Cu collisions at tJsnn = 
200 GeV, fgeo = 0.40 ± 0.047 [1(. Fitting ranges are limited to the bins with rii > 60, where rii is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


h 




^dyn 


x 2 

( n d) 


p- value 

[%] 


XLs/ n d with errors: 
quadrature statistical 
sum only 


0-5 


368271 


61.5 
±0.6 


19.79 
±0.01 


6.50 -lO -3 
±0.77 -10" 3 


1.13 
±0.02 


122.2 
4398.3 
(36) 





2.3 23.0 


5-10 


368869 


52.0 
±0.5 


16.74 
±0.01 


7.69 -10" 3 
±0.91 -10 -3 


1.13 
±0.02 


20.5 
613.9 

(30) 





0.39 3.9 


10-15 


369825 


42.3 
±0.4 


14.05 
±0.01 


9.46 -10~ 3 
±0.11 -10~ 2 


1.13 
±0.02 


16.2 
470.9 
(29) 





0.43 4.3 


15-20 


369964 


35.1 
±0.3 


11.77 
±0.01 


1.14 -10~ 2 
±0.13 -10 -2 


1.13 
±0.02 


11.4 
296.8 
(26) 





0.24 2.4 


20-25 


371752 


29.8 
±0.3 


9.80 
±0.01 


1.34 -10 -2 
±0.16 -10 -2 


1.13 
±0.02 


6.6 

370.4 

(23) 





0.20 2.0 


25-30 


368708 


25.6 
±0.3 


8.14 
±0.01 


1.56 -10" 2 
±0.18 -10 -2 


1.13 
±0.01 


42.7 
853.2 
(20) 





0.21 2.1 


30-35 


367869 


22.6 
±0.2 


6.72 
±0.005 


1.77 -10~ 2 
±0.21 -10 -2 


1.12 
±0.01 


126.4 
2274.4 





0.62 6.2 



(18) 



The results of this analysis show that this is too much. This is connected with the meaning of the formulation of a 
hypothesis. If the hypothesis is true, it means that the form of a distribution postulated by the hypothesis is exact 
for the whole population. Thus for the very large samples (as in all PHENIX cases) the measured distribution should 
be very close to that postulated. The performed analysis has shown that the PHENIX experimental errors are much 
bigger than the acceptable discrepancies (acceptable for these sample sizes). Therefore \ 2 an d Xls test statistics give 
the opposite answers. 
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TABLE VII: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Cu-Cu collisions at %/snn = 
62.4 GeV, f geo = 0.32 ± 0.063 0. Fitting ranges are limited to the bins with n; > 5, where rii is the number of events in the 
ith bin. 



Centrality 

[%] 


N 


k 


n 


1 / k>dyn 


bJdyn 


x 2 

i n d) 


p- value 

[%] 


xis/ n d with 
quadrature 
sum 


errors: 
statistical 
only 


0-5 


298182 


41.6 
±0.4 


13.35 
±0.01 


7.69 -lO -3 
±0.15 TO -2 


1.10 
±0.02 


9.3 
279.9 
(30) 





0.65 


6.5 


5-10 


307150 


26.5 
±0.2 


11.67 
±0.01 


1.21 TO -2 
±0.24 -10 -2 


1.14 
±0.03 


9.7 

290.7 
(30) 





0.78 


7.8 


10-15 


309874 


20.5 
±0.2 


9.90 
±0.01 


1.56 -10~ 2 

i n 01 i n — 2 

±0.31 • 10 


1.15 
±0.03 


9.3 
261.1 

(28) 





4.4 


44.0 


15-20 


312530 


17.8 
±0.1 


8.27 
±0.01 


1.80 TO -2 
±0.36 -10~ 2 


1.15 
±0.03 


26.0 
677.1 

(26) 





31.6 


316.1 


20-25 


312884 


16.0 
±0.1 


6.89 
±0.01 


1.99 -10~ 2 
±0.39 -10" 2 


1.14 
±0.03 


75.8 
1744.0 

(23) 





80.9 


809.3 



1.24 



CO 



dyn 



1.2- 



1.16 



1.12- 



1.08 



1.04 



O 



O 



<> 



<>A 



<> <> 



o 



A A 



<> j 
A <> 



A <> 



200 GeV Au-Au 
O PHENIX estimates 
A this work estimates 



CO 



~1 1 1 1 1 1 — 

100 200 300 400 



N 



part 



1.2 
dyn 

1.16- 
1.12- 
1.08 
1.04 



62.4 GeV Au-Au 
O PHENIX estimates 
A this work estimates 



<><> 



<> 



<> <> 



A A 



<> 



o 



<> 



<> 



<> 



~1 1 1 1 1 1 — 

100 200 300 400 



N 



part 



FIG. 7: Scaled variance for 200 GeV (left) and 62.4 GeV (right) Au-Au collisions. PHENIX estimates are from ilj. This work 
estimates are for the cases with ranges limited to the bins where m > 60, see Tables [Hi and II VI 
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TABLE VIII: Results of fitting multiplicity distributions measured by the PHENIX Collaboration in Cu-Cu collisions at 
75^7 = 62.4 GeV, f geo = 0.32 ± 0.063 Q. Fitting ranges are limited to the bins with rn > 60, where m is the number of 
events in the ith bin. 



Centrality 

[%] 


N 


k 


n 


1 / fcdyn 


bJdyn 


X 2 /nd 

x 2 

i n d) 


p- value 

[%] 


xls/ n d with 
quadrature 
sum 


errors: 
statistical 
only 


0-5 


298131 


42.0 
±0.5 


13.35 
±0.01 


7.62 -10~ 3 
±0.15 -10 -2 


1.10 
±0.02 


14.7 
411.9 

(28) 





0.67 


6.7 


5-10 


307061 


26.8 
±0.2 


11.66 
±0.01 


1.19 TO -2 
±0.24 -10 -2 


1.14 
±0.03 


19.7 
512.5 

(26) 





0.86 


8.6 


10-15 


309798 


20.7 

±U.2 


9.90 

±0.01 


1.54 -10" 2 

±0.30 ■ 10 


1.15 

i n no 
±0.03 


19.4 

4D0.0 

(24) 





0.38 


3.8 


15-20 


312434 


18.0 
±0.1 


8.27 
±0.01 


1.78 TO" 2 
±0.35 -10" 2 


1.15 
±0.03 


46.5 
976.4 
(21) 





0.40 


4.0 


20-25 


312758 


16.3 
±0.1 


6.89 
±0.01 


1.96 TO -2 
±0.39 -10" 2 


1.14 
±0.03 


118.1 
2243.4 
(19) 
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FIG. 8: Scaled variance for 200 GeV (left) and 62.4 GeV (right) Cu-Cu collisions. PHENIX estim ates are from Q. This work 
estimates are for the cases with ranges limited to the bins where rii > 60, see Tables IVTl and lyTTTl 
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Appendix A 

Dropping terms not depending on the parameters in Eq. (|17p , one obtains the following form for the log- likelihood 
function under consideration: 



lnL(Y \p,k) = p t X lnP(F i; p,fc) 

i=l 

Since the logarithm of the NBD is given by 



(Al) 



hiP(n;p, k) 

n 

= In (k+j — l) + nhi (1 — p) + k lnp — In (n!) 

i=i 

the necessary conditions for the existence of the maximum have the following form: 



(A2) 



ln£(Y | p, k) 



dp 

m 



Y- 



1 k 



1-p p 



= N 
= N 



1 m , rn 



= 



(A3) 



d_ 

Ok 



lni(Y | p, k) 



r Yi 



= N 



where the sum over j is if Yi = 0. 

From Eqs. (|A3I) and ([2]) one can obtain: 



N j2pr E 

=i L i=i 

m Yi 

E^E 



i 



i=l 



fc + j - 1 
1 

^ fc + j - 1 



hip 
lnp 



= 0, 



(A4) 



.„ , k(l-p) 
(N ch ) = -i ^ = n . 

P 

Expressing p as a function of k and (iV c h,) 

1 _ (Nch) | 1 ^ 
p A: 

and substituting it to Eq. (|A4[) the equation which determines fc is obtained: 



(A5) 



(A6) 



dk 



lni(Y|p, fc) 
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= N 



L i—1 j-1 J \ / 



= 



(A7) 



The above equation can be solved numerically. Having obtained k and substituting it into Eq. (|A6I) p is derived. 



Appendix B: Wilks's theorem 



Let X be a random variable with p.d.f f(X,9), which depends on parameters 6 = {9\, 
parameter space 9 is an open set in R d . For the set of N independent observations of X, X 
can defined the likelihood function 



?2,—,#d} G ©j where a 
= {X u X 2 , ...,X N }, one 



N 



Now consider iJg, a fc-dimensional subset of 0, A; < d. Then the maximum likelihood ratio can be defined as 



(Bl) 



max eggo L(X | 
maxg e e L(K \ ( 



(B2) 



This is a statistic because it does not depend on parameters 9 no more, in the numerator and the denominator there 
are likelihood function values at the ML estimators of parameters 6 with respect to sets H and 0, respectively. 

The Wilks's theorem says that under certain regularity conditions if the hypothesis Hq is true (i.e. it is true that 
6 e Hq), then the distribution of the statistic —2 In A converges to a \ 2 distribution with d — k degrees of freedom as 
TV — > oo 0,@. The proof can be found in Ref. Note that k = is possible, so one point in the parameter space 
(one value of the parameter) can be tested as well. 
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[12] Now Y x / Yj for i j and i,j = 1, 2, m. 



